Mon. Not. R. Astron. Soc. 000, [JHS] (2010) Printed 11 January 2013 (MN I^TeX style file v2.2) 



Transit timing variations in the HAT-P-13 planetary 
system 



o 



6 

> 

in 

in 
o 

(N 
O 



X 



Andras Pal^'^^f, Krisztian Sarneczky-*^, Gyula M. Szabo^'^, Attila Szing^, 
Laszlo L. Kiss^'^, Gyorgy Mezo^ and Zsolt Regaly^ 

^ Konkoly Observatory of the Hungarian Academy of Sciences, Konkoly Thege Miklos ut 15-17, Budapest, H-1121, Hungary, 
^ Department of Astronomy, Lordnd Eotvos University, Pdzmdny Peter setdny 1/A, Budapest H-llll, Hungary, 
^ Department of Experimental Physics and Astronomical Observatory, University of Szeged, 6720 Szeged, Hungary 
* Sydney Institute for Astronomy, School of Physics A28, University of Sydney, NSW 2006, Australia 



Accepted Received ... ; in original form 



ABSTRACT 

In this Letter we present observations of recent HAT-P-13b transits. The combined 
analysis of published and newly obtained transit epochs shows evidence for significant 
transit timing variations since the last publicly available ephemerides. Variation of 
transit timings result in a sudden switch of transit times. The detected full range of 
TTV spans ~ 0.015 days, which is significantly more than the known TTV events 
exhibited by hot Jupiters. If we have detected a periodic process, its period should be 
at least ~ 3 years because there arc no signs of variations in the previous observations. 
This argument makes unlikely that the measured TTV is due to perturbations by 
HAT-P-lSc. 
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1 INTRODUCTION 

Transit timing variations (TTVs) of transits are expected in 
transiting extrasolar planetary systems where the host star 
has more than one companion. This effect might reveal the 
presence of othe r companions with sm aller masses (see e.g. 
lAgol et al.ll2005l: ISteffen fc Agoll 120071 ') or co-orbital bodies 
as well ( Ford fc Holmanll2007l V On the other hand, if the sys- 
tem has two or more known planetary components, the mag- 
nitude of the TTV effects can be predi cted and/or used to 
characterize the sy stem more precisely (|Holman et al.|[2010l : 
ISteffen et all bOlOl 'l. Similarly, the lack o f TTVs can rule 
out o ther companions above a certain limit (|Csizmadia et all 
I2OIOI ). One should note that even two-body systems (i.e. the 
parent star and the planet) can show significant TTVs due 
to no n-gravitational physical processes (see e.g. iFabrvckvl 
|2008| ) and in the case of eccentric orbits, there are intrinsic 
tim ing variations due to the effects of general relativity (see 
e.g. iPal fc Kocsi'^l2008l ). In the case of two planetary com- 
panions, the magnitude of the TTV effects can be predicted 
analytically using the same methodology that is k nown for 
hierarchical stellar systems (|Borkovits et al.ll2003l ). and of 
course, by numerical integrations. 

At the time of its discovery, the planetary system or- 



biting the star HAT-P-13 was the only known multiple ex- 
trasolar planetary system that e xhibits at least on e transit- 
ing planet as well (HAT-P-13b. iBakos et al]|2009l ). Follow- 
up studies showed that the host star has a spin alignment 
with respect to th e orbital plane of the transiting planet 
IWinn et al.ll201(]| ') as well as this work also noticed a sig- 
nificant long-term drift in the radial velocity of the host 
star that might be interpret ed as a presence o f a third com- 
panion. Additional studies (|Szab6 et al.ll20l"ol ) reported the 
lack of primary transits of the long-period second compan- 
ion HAT-P-13c, and presented additional observations from 
the transits of the close-in plane t as well. As of thi s writing, 
no other (than lBakos et al.ll2009l : ISzab6 et al.ll2010l ) publicly 
available photometric data are available for this planetary 
system. Here we describe recent photometric measurements 
of the host star. The analysis of the light curves showed sig- 
nificant timing variations since the last public ephemerides. 
The structure of this Letter is as follows. In Section[2]we de- 
scribe details of the observations, data reduction and light 
curve modelling. The analysis of the results yielded by light 
curve fits is presented in Section [3] Section 3] discusses pos- 
sible scenarios for the observed TTVs and summarizes our 
results. 
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Figure 1. Photometry of the star HAT-P-13 on the nights of 
2010 December 27/28 (upper panel), December 30/31 (middle 
panel) and 2011 Januar 28/29 (lower panel), using the facilities 
of the Piszkesteto Mountain Station, Konkoly Observatory. Su- 
perimposed are the best-fit transit model light curves. Below the 
light curve plots, the fit residuals are also shown. Sec text for 
further details. 

2 OBSERVATIONS AND DATA REDUCTION 

We carried out photometric observations of the star HAT- 
P-13 with the telescopes of the Konkoly Observatory, lo- 
cated at the Piszkesteto Mountain Station on three sep- 
arate nights. On the night of 2010 December 27/28, we 
used the 60/90/180 cm Schmidt telescope, equipped with 
an Apogee ALTA-U 4k x 4k CCD camera. The following 
two events, on the night of December 30/31 and 2011 Jan- 



uar 28/29 has been partially observed (due to weather cir- 
cumstances on December 30 and due to dawn twilight on 
Januar 28) with the 100 cm Ritchey-Chretien-Coude (RCC) 
telescope equipped with a Princeton Instruments VersAr- 
ray 1.3k x 1.3k camera. These setups of optics and cam- 
eras yielded a square-shaped field-of-view with a size of 1.17 
and 0.11 degrees, for the Schmidt and RCC telescopes, re- 
spectively. At the first night we used Bessel I filter and an 
exposure time of 20 seconds (42 seconds cadence) while in 
the second and third night we used Cousins R filter and an 
exposure time of 60 and 25 seconds, respectively (here the 
readout time is relatively short, so this can be used as a 
cadence as well and on the last night the seeing was much 
better hence the shorter exposure time). 

In the following we describe the data reduction pro- 
cesses and some aspects of the light curve modelling. 



2.1 Data reduction 

Reduction of raw technical and scientific frames has been 
carried out with the pipeli nes built on the software pack- 
age described in EU (|2009h . The calibration procedure fol- 
lowed the standard way of dark subtraction and flat-field 
corrections using dark frames taken with similar exposures 
as the scientific images, (employing the tasks f icalib and 
f icombine). The source identification, astrometric solution 
and centroid coordinates were derived as follows. First, in- 
dividual point-like sources were extracted from each frame 
using the task f ista r. These lists of pro file coordinates were 
then cross- matched (jPal fc Bakosll2006l ) between the frames 
and a reference frame, yielding both a list of matched pairs 
and the spatial transformation between two images. The 
polynomial order of the spatial transformation were 3 and 
1 for the Schmidt and RCC fields, respectively (increasing 
these orders did not decrease the unbiased fit residuals). 
Then, using these initial differential astrometric solutions, 
the stellar pixel coordinates were refined using an indepen- 
dent profile centroid fit (done by the task fiphot) for the 
expected positions. These lists of refined profile coordinates 
were used to refine the spatial transformation as well. This 
is a rather relevant step because the point matching task 
(grmatch) sometimes misses a point and the exclusion of 
even a single point can yield a systematic deviation in the 
transformation coefficients whose error can propagate into 
the photometry errors (resulting in unexpected red noise). 
Following this two-stage astrometry, these refined coeffi- 
cients were used to transform the reference object list to 
the system of individual frames to perform aperture pho- 
tometry. Involving the task fiphot, aperture photometry 
was performed with annuli of appropriate sizes and a se- 
ries of apertures from which we picked the best one later 
on the modelling. Differential magnitudes were then com- 
puted from these raw instrumental magnitudes using sev- 
eral nearby comparison stars. The white-noise uncertain- 
ties of the individual photometric points were derive d from 
the p hoton noise, background noise and scintillation (jYound 
Il967h noise components of both the target star and compar- 
ison stars. 
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2.2 Light curve modelling 

Light curves obtained by the procedure (as described above) 
have been used to model the transit event and derive the re- 
spective parameters. In order to fit a model light curve to 
the observed d a ta, w e employed the analytic formulae of 
iMandel fc AgoJ (|2002l 'l. The effects of stellar limb darkening 
have been taken into account by fixing the quadratic limb 
dark ening coefficients u sing the stellar atmospheric parame- 
ters l|Bakos et al.ll2009l) . These c oefficients have been derived 
from the tables of lClaretl (|200(]| '). Since many of the regres- 
sion methods used by us require the knowledge of paramet- 
ric derivati ves o f the m odel function, we used the formulae 
provided bv lPall (|2008l ) to compute these partial derivatives. 

The light curve models have been parameterised using 
p = Rp/Ri,, the planet-to-star size ratio, b^, the square 
of the impact parameter, the quantity C^/R*, related to 
the duration of the transit as C/^* = 2/(r3.5 — T1.5) and 
Tc = (Ti.s + r3.5)/2, the time instance at the center of 
the transit (here Ti.s and T3.5 denotes the time instances 
when the center of the planet disc intersects the lim b of 
the st ar inwards and outwards, respectively). See also I Pall 
l|2008h for further details on the advantages of this param- 
eterisation. In addition to these parameters, we fitted si- 
multaneously the out-of-transit magnitude, a linear trend 
in the out-of-transit magnitude (correcting to the effects of 
changing airmasses) and a linear decorrelation parameter 
against the variations resulted by the changes in the profile 
FWHM. In the case of the observation on 2010 December 
30/31, we kept 7?p/_R*, and C/R* fixed, since this light 
curve does not cover the entire transit event. We emp loyed 
the M arkov Chain Monte Carlo method (MCMC; see iFordI 
I2OO4I ) to obtain the a posteriori distribution of the fit pa- 
rameters from which the best fit values, uncertainties and 
correlations can easily be derived. The best fit values for 
the light curve shape parameters for the event 2010 Decem- 
ber 27/28 are Rp/R^, = 0.0880 ± 0.0031, 6^ = 0.544 ± 0.086 
and C/R* = 17.07 ± 0.28. These values agree s well with the 
previously published ones (jBakos et al ] |2009l ). Fig. [U shows 
the light curves with these best fit model functions superim- 
posed while the best fit values for the transit center instances 
are listed in Table [1] This table also lists the previously 
published individu al tr ansit event instanc es as well, from 
iBakos et all (|2009l ) and ISzabo et al.l (|2010l ). The fit residu- 
als has been compared to the expected residual value that 
can be computed from the photon noise, background noise 
and scintillation of the star. We find that the fit residuals 
for the observation of 2010 December 27/28 were only larger 
by a factor of 1.05, that indicates a very good quality light 
curve with small systematic variations (red noise). Indeed, 
even the raw instrumental magnitudes (i.e. the magnitude 
of the target star without subtracting any comparison star) 
varied in a range of only 0.02 mag. However, the analysis 
of the observation from 2010 December 30/31 yielded a fit 
residual that is larger by a factor of « 2 than the expected 
and as it clear from the lower plot on Fig. [1] the residu- 
als contain significant amount of red noise. This was due 
to both weather circumstances and significant flat fleld im- 
age structures in the vicinity of HAT-P-13. Therefore we 
conservatively scaled up the derived errors by that factor 
of 2 (and this is the value that has been reported in Ta- 
ble [1]). The observation carried out on 2011 January 28/29 



Table 1. Transit center instances from previous works and re- 
cent measurements. Data for the events —68 . . . -1-62 are available 
from Bakos et al. (2009), for the events 124 and 161 data are 
taken from Szabo et al. (2010). The last three fields are the new 
measurements presented in this Letter. The event enumeration is 
from the same reference epoch as in Bakos et al. (2009). See text 
for further details. 



Event 


BJD 


Error 


-68 


2454581.62406 


0.00122 


-1 


2454777.01287 


0.00100 





2454779.92953 


0.00063 


1 


2454782.84357 


0.00155 


24 


2454849.92062 


0.00075 


35 


2454882.00041 


0.00150 


62 


2454960.73968 


0.00178 


124 


2455141.55220 


0.00100 


161 


2455249.45080 


0.00200 


267 


2455558.56265 


0.00098 


268 


2455561.48379 


0.00400 


278 


2455590.64486 


0.00179 



was partial: the out-of-transit part of the light curve after 
the transit was not measured, however, the observed egress 
can be incorporated into the light curve analysis. This fit has 
been done similarly as for the night of December 27/28 and 
yielded Rp/R^ = 0.0820 ± 0.0078, b^ = 0.471 ± 0.214 and 
C,/Ri, = 16.88 ± 0.48. These values also agrees well (within 
uncertainties) both with the fit of December 27/28 data and 
the published values. For the center of the transit, we obtain 
the value found in Table [T] For this analysis, the uncertain- 
ties yielded by the MCMC fit have been multiplied by 1.5 
in order to take into account the red noise components: this 
value is the ratio of the fit residuals and the expected for- 
mal photometric errors (those have been derived from the 
photon and background noise and the scintillation) . The un- 
certainties reported above for the geometric parameters and 
listed in Table [T] are these increased values. 



3 TIMING VARIATIONS 

Using the data of iBakos et all (|2009l ') and [S zabo et al.l 
(2010| ), the subsequent observed transit events can be mod- 
elled easily with a strictly periodic assumption: the formal 
errors reported in these papers yield an unbiased resid- 
ual of xo = 9.2 for these 7 degrees of freedom. By defin- 
ing the event A'tr = as the reference epoch, this lin- 
ear model yields E = 2454779.92978 ± 0.00049 (BJD) and 
P = 2.9162953±0.0000085 days while Corr(S,P) = -0.406. 
These ephemerides gives an uncertainty of At = 0.0021 days 
for the predictions of the discussed events (A^tr = 267, 268 
and 278 for the nights 2010 December 27/28, 30/31 and 
January 28/29). As it is clear from Fig. [2l these three re- 
cent measurements outlie from the linear fit by 8.4-cr, 3.3-(t 
and 5.5-(T, respectively. We note here that the significance 
of the first (December 27/28) measurement is determined 
mostly by the uncertainty of the predictions while the sig- 
nificance of the second one is determined by the uncertainty 
of that particular measurement. All in all, both values can 
be treated as a significant deviation. 

Due to the significance of the results, i.e. the deviance of 
the transit center instances from the predictions, we applied 
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Figure 2. Transit timing variations as the function of the event 
number (cycles since event at the night of 2008 November 9/10). 
The points with the errorbars represent the individual measure- 
ments and the respective uncertainties as listed in Table. [T] The 
thick solid line is the best fit linear ephemerides to the first 9 
data points while the thinner dashed lines above and below it 
represents the l-cr errorbars of the best linear fit. 



several independent analysis methods. These methods in- 
clude the minimisation aleorithm based on the downhill sim- 



plex method (jPress et al.lll992l 'l. as well as a one-parameter 
scan. In the latter test, the light curve shape parameters 
{Rp/Ri,, and C/R*) have been fixed and only the Tc has 
been varied while the out-of-transit magnitude and other 
decorrelation parameters (see above) have been minimised 
using the classic linear least squares method. Additionally, 
we compared the statistical covariances with the covariances 
yielded by linear error propagation analysis. All of these 
methods have confirmed the results and uncertainties of the 
Markov Chain Monte Carlo analysis presented in the previ- 
ous section. In addition to these tests, we checked other data 
acquistion logs to exclude other sources of possible system- 
atic effects. The filesystem metadata information has been 
compared with the FITS keywords as well as the logs of the 
NTP daemon were also checked (this showed an approxi- 
mately 0.002 seconds of standard deviation from the time 
synchronization servers) . 



4 DISCUSSION 

There are a few examples for extrasolar planets proven to 
exhibit TTV. WASP-3b has been found to display a full 
TTV amplitude of 0.004 day period. As an explanation, a 
third planet has been inferred on a resona nce orbit that may 
be eit her an inner or outer perturber (|Macieiewski et al] 
l2010al ). The TTV of WASP-5 is about 0.0023 day, and a 
set of different perturber solutions were invok ed that can 
equiv alent ly well explain the observed TTV (|Fukui et al.l 
I2OIOI ). WASP-10 exhibits a similar TTV to that of WASP- 
3 and WASP-5, its TTV has a full amplitude of 0.002 day 
with 5.47 day period. The perturber has been suggested to 
have 5Mj mass and orbit in 5 : 3 mean m otion resonance 
with WASP-lOb (|Macieiewski et all [2010bl ') . Other sources 
than perturbing planets are also known as sources of TTV , 
such as exomoons orbiting the planets (jSzabo et al.ll200^ : 



ISimon. Szatmarv fc SzabdIbOOTi ). However, the transit tim- 
ing variations are in the order of a few seconds/minutes in 
the case of plausible moons, and the measured TTV of HAT- 
P-13b is not compatible with this predicted TTV amplitude. 

The detected O - C of HAT-P-13b is in the order of 
0.01 days, significantly more than that of WASP-3, WASP- 
5 or WASP-10. Kepler-9b,c represent an example for such 
large TTV amplit ude exhibited by Ne ptunes that orbit in 
strong resonance (jHolman et al.ll2010l l. However, HAT-P- 
13b is not very similar to the case of Kepler-9b,c. The plan- 
ets in the Kepler-9 system exhibit a continuous variation of 
TTV with a period of few hundred days, while HAT-P-13b 
seemed to lack TTV during the first three years of follow-up 
observations, and then a sudden switch of TTV is observed. 

This observed behaviour and the TTV amplitude of 
HAT-P-13b may be explained by perturbations by long 
period planets on no n-resonant orbits. For instance, in 
iBorkovits et al.l (|2010l ). a model planet similar to CoRoT- 
9 exhibits a large switch in O — C of the inner planet in 
less than lyear. Assuming an outer perturber of 5 Afj mass, 
10, 000 day orbital period and e = 0.7 eccentricity, the pre- 
dicted amplitude results to be 0.013 days (see Fig. 6 in the 
cited paper). This value is quite similar to the O — C varia- 
tion that we report for HAT-P-13, suggesting that perturba- 
tions on the time-scale of the orbital period of a long-period 
perturber can explain our inferred TTV of HAT-P-13b. The 
detailed modeling is beyond the information content of the 
currently available dataset: there would be at least 6 orbital 
parameters to explain with the amplitude and the timing of 
the switch of the O — C diagram, leading to a highly un- 
constrained problem. It may be unlikely that we observed 
perturbations of HAT-P-13c, because this non-resonant per- 
turber would lead to TTV variations with 428.5 day ampli- 
tude and additional secular components. This period is cov- 
ered by the current dataset and there are no signs for this 
period in the deduced TTV. Eventually, the most urgent 
task is deriving the period of the TTV shown by HAT-P- 
13b before predictions can be given to the nature of the 
perturber. 
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